args <- commandArgs(TRUE)
if(length(args) != 6){
    cat ("Usage:\n\tRscript quality.R <Base_quality_value_distribution_by_read_position_1.txt> <Base_quality_value_distribution_by_read_position_2.txt> <output raw png file> <output clean png file> <read1 trim> <read2 trim>\n")
    q(save="no")
}
line1 <- as.numeric(system(paste("wc -l ",args[[1]],"|awk '{print $1}'"),intern=T))
line1 <- (line1+as.numeric(args[[5]])-4)/2
rawrt1 <- read.table(args[[1]],skip=2,nrows=line1)
cleanrt1 <- read.table(args[[1]],skip=line1+4)
line2 <- as.numeric(system(paste("wc -l ",args[[2]],"|awk '{print $1}'"),intern=T))
line2 <- (line2+as.numeric(args[[6]])-4)/2
rawrt2 <- read.table(args[[2]],skip=2,nrows=line2)
cleanrt2 <- read.table(args[[2]],skip=line2+4)
rawrt1 <- rev(rawrt1)
cleanrt1 <- rev(cleanrt1)
rawrt2 <- rev(rawrt2)
cleanrt2 <- rev(cleanrt2)
pos <- c(rawrt1$V1,rawrt2$V1+max(rawrt1$V1))
rawmean <- c(rawrt1[[6]],rawrt2[[6]])
rawmedian <- c(rawrt1[[5]],rawrt2[[5]])
rawlower <- c(rawrt1[[4]],rawrt2[[4]])
rawupper <- c(rawrt1[[3]],rawrt2[[3]])
raw10th <- c(rawrt1[[2]],rawrt2[[2]])
raw90th <- c(rawrt1[[1]],rawrt2[[1]])
cleanmean <- c(cleanrt1[[6]],cleanrt2[[6]])
cleanmedian <- c(cleanrt1[[5]],cleanrt2[[5]])
cleanlower <- c(cleanrt1[[4]],cleanrt2[[4]])
cleanupper <- c(cleanrt1[[3]],cleanrt2[[3]])
clean10th <- c(cleanrt1[[2]],cleanrt2[[2]])
clean90th <- c(cleanrt1[[1]],cleanrt2[[1]])
rawstat <- matrix(c(raw10th,rawlower,rawmedian,rawupper,raw90th),line1+line2,5)
rawstat <- t(rawstat)
raw <- list(stats=rawstat,n=rep(Inf,line1+line2))
cleanstat <- matrix(c(clean10th,cleanlower,cleanmedian,cleanupper,clean90th),line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]),5)
cleanstat <- t(cleanstat)
clean <- list(stats=cleanstat,n=rep(Inf,line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]])))
png(file=args[[3]],height=600,width=1500,antialias='default')
bxp(raw,main="Base quality along reads",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(0,42))
lines(rawmean,col="red",lwd=2)
abline(v=line1+.5,col="blue",lty=4,lwd=2)
axis(side=1, seq(0,line1+line2,20), tcl=0.2, labels=FALSE)
axis(side=2, seq(0,40,10), tcl=0.2, labels=FALSE)
axis(side=3, seq(0,line1+line2,20), tcl=0.2, labels=FALSE)
axis(side=4, seq(0,40,10), tcl=0.2, labels=FALSE)
mtext(seq(0,line1+line2,20),side=1,las=1,at=seq(0,line1+line2,20),line=0.3,cex=1.5)
mtext(seq(0,40,10),side=2,las=1,at=seq(0,40,10),line=0.3,cex=1.5)
mtext("Position along reads",side=1, line=2, at=(line1+line2)/2, cex=1.5)
mtext("Quality",side=2, line=2.5, at=21, cex=1.5)
d <- dev.off()
png(file=args[[4]],height=600,width=1500,antialias='default')
bxp(clean,main="Base quality along reads",xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(0,42))
lines(cleanmean,col="red",lwd=2)
abline(v=line1-as.numeric(args[[5]])+.5,col="blue",lty=4,lwd=2)
axis(side=1, seq(0,line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]),20), tcl=0.2, labels=FALSE)
axis(side=2, seq(0,40,10), tcl=0.2, labels=FALSE)
axis(side=3, seq(0,line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]),20), tcl=0.2, labels=FALSE)
axis(side=4, seq(0,40,10), tcl=0.2, labels=FALSE)
mtext(seq(0,line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]),20),side=1,las=1,at=seq(0,line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]),20),line=0.3,cex=1.5)
mtext(seq(0,40,10),side=2,las=1,at=seq(0,40,10),line=0.3,cex=1.5)
mtext("Position along reads",side=1, line=2, at=(line1+line2-as.numeric(args[[5]])-as.numeric(args[[6]]))/2, cex=1.5)
mtext("Quality",side=2, line=2.5, at=21, cex=1.5)
d <- dev.off()

